\begin{section}{Apéndices}
\newcommand{\real}{\ensuremath{\mathbb{R}}}

\begin{subsection}{Apéndice A - Enunciado} % Enunciado
\begin{centering}
\bf Laboratorio de M\'etodos Num\'ericos - Primer cuatrimestre de 2011 \\
Trabajo Pr\'actico N\'umero 3: CAD - M\'as que splines \\
\end{centering}

\vskip 25pt
\hrule
\vskip 11pt

% {\bf Enunciado}

Los programas de dise\~no asistido por computadora (CAD) son herramientas fundamentales para ingenieros, arquitectos, dise\~nadores, artistas y animadores. Sus interfaces gr\'aficas esconden un sinn\'umero de complicadas operaciones. Un ejemplo b\'asico de tales operaciones es la tarea de seleccionar un punto cualquiera de una curva y moverlo a una nueva posici\'on deformando la curva. El punto seleccionado se ingresa usualmente mediante un dispositivo apuntador (\emph{mouse} o tableta digitalizadora) interactuando con la interfaz, y en general el punto ingresado est\'a meramente \emph{cerca} de la curva.

En este trabajo pr\'actico se deber\'an dise\~nar algoritmos e implementar un programa que, dadas las coordenadas $(\bar{x}_i,\bar{y}_i) \in \real^2$ de una serie de puntos de control ($i = 1..n$), construya una spline natural\footnote{Esto es suficiente para nuestro TP, pero en realidad los sistemas de CAD utilizan m\'as frecuentemente otros mecanismos para obtener, describir y manipular curvas y superficies.} param\'etrica que pase por los puntos en el orden dado. Adem\'as, dadas las coordenadas $(x^*,y^*) \in \real^2$ de un punto cercano a la curva, calcule el punto de la curva m\'as pr\'oximo y construya una nueva spline natural resultante de modificar la spline original de forma que ahora adem\'as pase por la nueva posici\'on $(\bar{x}^*,\bar{y}^*) \in \real^2$ del punto seleccionado. 
El procedimiento descripto se muestra en la figura con la spline original dibujada en l\'inea de trazos y la spline deformada en l\'inea continua.
\begin{figure}[H] \centering
	\subfloat{\includegraphics[scale=0.06]{./graficos/img_enunciado.eps}}
	%~ \caption{Spline}
\end{figure}
El programa deber\'a trabajar las splines como curvas param\'etricas en $\real^2$. Dadas las coordenadas de los puntos de control existen varias estrategias para definir la parametrizaci\'on. Algunas de las parametrizaciones com\'unmente utilizadas son:
%~ \vspace{-15pt}
\begin{description}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
 \item[Uniforme:] la variaci\'on del par\'ametro es igual entre cualquier par de puntos de control consecutivos;
 \item[\emph{Chord-length}:] la variaci\'on del par\'ametro entre dos puntos de control consecutivos es proporcional a la distancia entre los mismos;
 \item[Centr\'ipeta:] la variaci\'on del par\'ametro es proporcional a la ra\'iz cuadrada de la distancia entre los puntos de control\footnote{Este m\'etodo fue propuesto por Eugene Lee en \emph{Choosing nodes in parametric curve interpolation}, Computer-Aided Design 21, 1989.}.
\end{description}
%~ \vspace{-15pt}
En este trabajo pr\'actico deber\'an utilizar alguna de estas parametrizaciones. Opcionalmente podr\'an implementar las restantes y comparar los resultados obtenidos con las tres variantes.

Adem\'as, el programa deber\'a conservar el valor del par\'ametro que le corresponde a cada punto de control y al punto seleccionado, antes y despu\'es de moverlo. 
\newpage
{\bf Preguntas:}
%~ \vspace{-15pt}
\begin{enumerate}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
 \item >Depende la forma de la curva de la elecci\'on de la parametrizaci\'on? 
 \item >Cambia la forma de la curva si en lugar de deformar la curva conservando la parametrizaci\'on el programa la recalcula al mover el punto?
 \item (Opcional) >C\'omo cambia la forma de la curva seg\'un la condici\'on de borde usada (natural, sujeto, \textit{not-a-knot}, etc.)? 
 \item (Opcional) Si se quiere redibujar cont\'inuamente la curva mientras el usuario mueve el punto seleccionado, >c\'omo se puede calcular esto m\'as eficientemente?
 \item (Opcional) Luego de mover el punto seleccionado, >cambia toda la curva (\emph{control global}) o solamente una parte (\emph{control local})? >Qu\'e consecuencias puede tener esto?
 \item (Opcional) Si el intervalo del par\'ametro se muestrea uniformemente, >los puntos resultantes quedan espaciados uniformemente? >Qu\'e otras alternativas de muestreo ser\'ian apropiadas? 
 \item (Opcional) Si se necesitara que las longitudes de curva entre puntos consecutivos fueran todas iguales, >c\'omo deber\'ia muestrearse?
\end{enumerate}

{\bf Archivos de entrada / salida}

La entrada de datos se realizar\'a mediante un archivo de texto con el siguiente formato:
%~ \vspace{-15pt}
\begin{itemize}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
\item En la primera l\'inea figurar\'a el n\'umero $n$ de puntos de control utilizados para definir la spline y, separado por espacio, el n\'umero $m$ de puntos de muestreo de la spline. 
\item En las siguientes $n$ l\'ineas figurar\'an las coordenadas $\bar{x}$ e $\bar{y}$ de cada punto de control separadas por espacio.
\item Una l\'inea en blanco
\item Una l\'inea con las coordenadas $x^*$ e $y^*$ del punto pr\'oximo a la curva, separadas por espacio.
\item Una l\'inea en blanco
\item Una l\'inea con las coordenadas $\bar{x}^*$ e $\bar{y}^*$ de la nueva posici\'on del punto, separadas por espacio.
\end{itemize}

La salida de datos estar\'a dada por un archivo de texto con el siguiente formato:
%~ \vspace{-15pt}
\begin{itemize}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
\item En la primera l\'inea figurar\'a el n\'umero $m$ de puntos muestreados. 
\item En las siguientes $m$ l\'ineas figurar\'an las coordenadas $x$ e $y$ de cada punto muestreado en la spline original, separadas por espacio. Estos puntos corresponder\'an a un muestreo uniforme del rango del par\'ametro e incluir\'an los extremos. De esta forma, probablemente este conjunto de puntos no incluya los puntos de control originales.
\item Una l\'inea en blanco
\item Una l\'inea con las coordenadas del punto en la curva original m\'as pr\'oximo al punto ingresado, separadas por espacio.
\item Una l\'inea en blanco
\item En las siguientes $m$ l\'ineas figurar\'an las coordenadas $x$ e $y$ de cada punto muestreado en la spline deformada, separadas por espacio.
\end{itemize}


\vskip 30pt
\hrule
\vskip 11pt

{\bf Entregas parciales}
%~ \vspace{-15pt}
\begin{description}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
 \item[10 de junio:] Implementaci\'on del c\'alculo de splines y funciones de entrada/salida.
 \item[17 de junio:] Implementaci\'on completa y verificaci\'on.
\end{description}

{\bf Entrega Final}
%~ \vspace{-15pt}
\begin{description}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
 \item[Formato Electr\'onico:] 23 de junio de 2011, hasta las 23:59 hs, a la direcci\'on: 

  {\emph{metnum.lab2011@gmail.com}}
 \item[Formato f\'isico:] 24 de junio de 2011, de 17 a 21 hs.
\end{description}
\end{subsection}

\newpage
\begin{subsection}{Apéndice B - Código} % Código fuente relevante y algunos resultados específicos
%~ \newpage
\begin{subsubsection}{main}
\begin{lstlisting}
#include "param.h"
#include "spline.h"
#include "evaluar.h"

void cbool(bool a)
{
	if (a) 
		{cout << "True";}
	else
		{cout << "False";}
	
	return;
}

/* Variables Globales. */
long int cant_puntos 	= 1000;
bool logs 		= true;
bool mostrar_splines 	= false;
bool mostrar_puntos 	= false;
bool mostrar_evaluacion	= false;
bool mostrar_graficos 	= true;

/* ./main entrada.in salida.out nombre_fig */
int main(int argc, char** argv)
{	
	string nombre_entrada, nombre_salida, nombre_fig;
	
	if (argc == 4)
	{
		nombre_entrada = argv[1];		
		nombre_salida = argv[2];		
		nombre_fig = argv[3];		
	}
	else
	{
		if (argc == 3)
		{
			nombre_entrada = argv[1];		
			nombre_salida = argv[2];		
			nombre_fig = argv[1];		
		}
		else
		{
			if (argc == 2)
			{
				string nombre = argv[1];
				nombre_entrada = nombre + ".in";
				nombre_salida = nombre + ".out";
				nombre_fig = nombre;
			}	
			else 
			{
				nombre_entrada = "spline.in";
				nombre_salida = "spline.out";		
				nombre_fig = "spline";		
			}
		}
	}
	
	ifstream entrada(nombre_entrada.c_str());
	ofstream salida(nombre_salida.c_str());	
	
	cout.precision(20);
	
	int n, m;
	/* n: # puntos de control 
	 * m: # puntos de muestreo. */
	
	entrada >> n >> m;
	
	if (logs)
		cout << "n : " << n << endl;
		
	/* Hacer que la cantidad de puntos en los que se va a 
	 * dividir la curva dependa de n. */
	cant_puntos = cant_puntos * n;
	
	if (n == 1)
	{
		cerr << "Debe ingresar mas de un punto." << endl;
		exit(1);
	}
		
	/* Coordenadas (x,y) de los puntos de control. */
	vPunto coordenadas;	
	forn(i, n)
	{
		double X, Y;
		entrada >> X >> Y;
		coordenadas.push_back(Punto(X, Y));		
	}
	
	if (logs)
	{
		cout << endl << "Coordenadas (x,y) de los puntos" << endl;
		
		forn(i, n-1)
			cout << "(" << coordenadas[i].x << ", " << coordenadas[i].y << "), ";
		cout << "(" << coordenadas[n-1].x << ", " << coordenadas[n-1].y << ")" << endl;
	}
	
	double X, Y;
	entrada >> X >> Y;
	Punto click(X, Y);				
	entrada >> X >> Y;
	Punto nuevo(X, Y);
	
	if (logs)
	{
		cout << endl << "Click: (" << click.x << ", " << click.y << ") Nueva posicion: (" << nuevo.x << ", " << nuevo.y << ")" << endl;
	}
		
	/* Fin lectura de archivo de entrada. */
	
	vdouble param (n);
	double  var = 1.0;
	//~ param_Uniforme(param, n, var);
	param_Chord(coordenadas, param, n, var);
	//~ param_Centripeta(coordenadas, param, n, var);
	enRango(param, n);
	
	if (logs)
	{
		cout << endl << "Parametrizacion " << endl;	
		forn(i, n)
			cout << param[i] << " ";	
		cout << endl;
	}
	
	/* Armar nuevos puntos de la forma (s_i, x(s_i)) y (s_i, y(x_i)). */
	vPunto puntos_x (n);
	vPunto puntos_y (n);
	param_puntos(coordenadas, param, puntos_x, puntos_y, n);
	
	if (mostrar_puntos)
	{
		cout << endl << "Puntos en x" << endl;
		forn(i, n-1)
			cout << "(" << puntos_x[i].x << ", " << puntos_x[i].y << "), ";
		cout << "(" << puntos_x[n-1].x << ", " << puntos_x[n-1].y << ")" << endl;
		
		cout << endl << "Puntos en y" << endl;
		forn(i, n-1)
			cout << "(" << puntos_y[i].x << ", " << puntos_y[i].y << "), ";
		cout << "(" << puntos_y[n-1].x << ", " << puntos_y[n-1].y << ")" << endl;
	}
	
	/* Armar los splines. */
	vSpline spline_x (n-1);
	vSpline spline_y (n-1);
	
	/* Cuando n es 2 se resuelve por separado. */
	if (n == 2)
	{
		spline_x[0].a = puntos_x[0].y;
		spline_x[0].b = (puntos_x[1].y - puntos_x[0].y)/(puntos_x[1].x - puntos_x[0].x);
		spline_x[0].c = 0;
		spline_x[0].d = 0;
		spline_x[0].x_i = puntos_x[0].x;
		spline_x[0].min = puntos_x[0].x;
		spline_x[0].max = puntos_x[1].x;
		
		spline_y[0].a = puntos_y[0].y;
		spline_y[0].b = (puntos_y[1].y - puntos_y[0].y)/(puntos_y[1].x - puntos_y[0].x);
		spline_y[0].c = 0;
		spline_y[0].d = 0;
		spline_y[0].x_i = puntos_y[0].x;
		spline_y[0].min = puntos_y[0].x;
		spline_y[0].max = puntos_y[1].x;		
	}
	else
	{
		interpolar(puntos_x, spline_x, n);
		interpolar(puntos_y, spline_y, n);	
	}
	
	if (mostrar_splines)
	{
		cout << endl << "Spline en x" << endl;	
		forn(i, n-1)
		{		
			double punto = spline_x[i].x_i;
			cout << "S_" << i << "(x) = " << spline_x[i].d << "(x-" << punto << ")^3 + ";
			cout << spline_x[i].c << "(x-" << punto << ")^2 + ";
			cout << spline_x[i].b << "(x-" << punto << ") + ";
			cout << spline_x[i].a;
			cout << " x en [" << spline_x[i].min << ", " << spline_x[i].max << "]" << endl;
		}
		
		cout << endl << "Spline en y" << endl;		
		forn(i, n-1)
		{		
			double punto = spline_y[i].x_i;
			cout << "S_" << i << "(y) = " << spline_y[i].d << "(y-" << punto << ")^3 + ";
			cout << spline_y[i].c << "(y-" << punto << ")^2 + ";
			cout << spline_y[i].b << "(y-" << punto << ") + ";
			cout << spline_y[i].a;
			cout << " y en [" << spline_y[i].min << ", " << spline_y[i].max << "]" << endl;
		}		
	}
	
	if (logs)
	{
		bool res;
		res = verificar(puntos_x, spline_x, n);
		if (!res)
			cerr << "Hubo un problema al calcular el spline en x" << endl;	
		res = verificar(puntos_y, spline_y, n);
		if (!res)
			cerr << "Hubo un problema al calcular el spline en y" << endl;
	}
	
	/* Calcular uniformemente los puntos de muestreo. */
	vdouble muestreo (m);
	muestreo_Uniforme(muestreo, m);
	
	if (logs)
	{
		cout << endl << "Puntos de muestreo" << endl;
		forn(i, m)
			cout << muestreo[i] << " ";
		cout << endl;
	}
	
	/* Obtener coordenadas (x,y) de cada punto de muestreo segun la spline. */
	vPunto muestreados (m);
	evaluar(muestreo, muestreados, spline_x, spline_y, n, m);
	
	if (mostrar_evaluacion)
	{
		cout << endl << "Muestreados" << endl;
		forn(i, m)	
			cout << muestreados[i].x << " " << muestreados[i].y << endl;
	}
	
	if (mostrar_graficos)
	{
		graficar(cant_puntos, nombre_fig, nombre_fig + " curva inicial.", spline_x, spline_y, n);	
	}
	
	/* Calcular el punto perteneciente al spline mas proximo al punto 
	 * en el cual el usuario hizo click. */
	Punto proximo;
	double t;
	t = mas_Proximo(cant_puntos, spline_x, spline_y, click, proximo, n);
	
	if (logs)
	{
		cout << endl << "Punto mas proximo" << endl;
		cout << proximo.x << " " << proximo.y << " t: " << t << " Distancia al click: " << fabs(distancia(click, spline_x, spline_y, t, n)) << endl;
	}
	
	/* Nuevos puntos para x y para y. */
	Punto nuevo_x;
	Punto nuevo_y;
	vPunto nuevos_puntos_x;
	vPunto nuevos_puntos_y;
	
	nuevo_x.x = t;
	nuevo_x.y = nuevo.x;
	nuevo_y.x = t;
	nuevo_y.y = nuevo.y;
	
	int pos = ubicacion(nuevo_x, puntos_x);
	agregar_punto(nuevo_x, pos, puntos_x, nuevos_puntos_x);
	agregar_punto(nuevo_y, pos, puntos_y, nuevos_puntos_y);
	
	if (mostrar_puntos)
	{
		cout << endl << "Nuevos puntos en x" << endl;
		forn(i, n)
			cout << "(" << nuevos_puntos_x[i].x << ", " << nuevos_puntos_x[i].y << "), ";
		cout << "(" << nuevos_puntos_x[n].x << ", " << nuevos_puntos_x[n].y << ")" << endl;
		
		cout << endl << "Nuevos puntos en y" << endl;
		forn(i, n)
			cout << "(" << nuevos_puntos_y[i].x << ", " << nuevos_puntos_y[i].y << "), ";
		cout << "(" << nuevos_puntos_y[n].x << ", " << nuevos_puntos_y[n].y << ")" << endl;
	}
		
	/* Agregar una posicion a los splines. */	
	spline_x.push_back(Spline());
	spline_y.push_back(Spline());
	
	/* Interpolar nuevamente usando los nuevos puntos. */
	interpolar(nuevos_puntos_x, spline_x, n+1);
	interpolar(nuevos_puntos_y, spline_y, n+1);	
	
	if (mostrar_splines)
	{
		cout << endl << "Spline en x" << endl;	
		forn(i, n)
		{		
			double punto = spline_x[i].x_i;
			cout << "S_" << i << "(x) = " << spline_x[i].d << "(x-" << punto << ")^3 + ";
			cout << spline_x[i].c << "(x-" << punto << ")^2 + ";
			cout << spline_x[i].b << "(x-" << punto << ") + ";
			cout << spline_x[i].a;
			cout << " x en [" << spline_x[i].min << ", " << spline_x[i].max << "]" << endl;
		}
		
		cout << endl << "Spline en y" << endl;		
		forn(i, n)
		{		
			double punto = spline_y[i].x_i;
			cout << "S_" << i << "(y) = " << spline_y[i].d << "(y-" << punto << ")^3 + ";
			cout << spline_y[i].c << "(y-" << punto << ")^2 + ";
			cout << spline_y[i].b << "(y-" << punto << ") + ";
			cout << spline_y[i].a;
			cout << " y en [" << spline_y[i].min << ", " << spline_y[i].max << "]" << endl;
		}
	}
	
	if (logs)
	{
		bool res;
		res = verificar(nuevos_puntos_x, spline_x, n+1);
		if (!res)
			cerr << "Hubo un problema al calcular el spline en x" << endl;	
		res = verificar(nuevos_puntos_y, spline_y, n+1);
		if (!res)
			cerr << "Hubo un problema al calcular el spline en y" << endl;
	}
			
	vPunto deformados (m);
	evaluar(muestreo, deformados, spline_x, spline_y, n+1, m);	
	
	if (mostrar_evaluacion)
	{
		cout << endl << "Deformados" << endl;
		forn(i, m)	
			cout << deformados[i].x << " " << deformados[i].y << endl;
	}
		
	/* Escritura archivo de salida. */
	salida << m << endl;
		
	forn(i, m)	
		salida << muestreados[i].x << " " << muestreados[i].y << endl;	
	
	salida << endl;	
	salida << proximo.x << " " << proximo.y << endl;
	salida << endl;
	
	forn(i, m)	
		salida << deformados[i].x << " " << deformados[i].y << endl;
	
	return 0;
}

\end{lstlisting}
\end{subsubsection}

%~ \begin{subsubsection}{main}
%~ \begin{lstlisting}
%~ 
%~ \end{lstlisting}
%~ \end{subsubsection}

\begin{subsubsection}{param}
\begin{lstlisting}
#include "param.h"

double dist(Punto p1, Punto p2)
{
	double X = (p1.x + p2.x);
	double Y = (p1.y + p2.y);
	
	double res = sqrt(X*X + Y*Y);
	
	return res;
}

/* Parametrizaciones. */
void param_Uniforme(vdouble& param, int n, double variacion)
{
	double valor = 0;
	
	forn(i, n)
	{
		param[i] = valor;
		valor += variacion;
		cout << param[i] << " ";
	}
	
	return;
}

void param_Chord(vPunto& coordenadas, vdouble& param, int n, double variacion)
{
	double valor = 0;
	
	forn(i, n-1)
	{
		param[i] = valor;
		valor += variacion * dist(coordenadas[i], coordenadas[i+1]);
	}
	
	/* n siempre es mayor a 0 */
	param[n-1] = valor;
	
	return;
}

void param_Centripeta(vPunto& coordenadas, vdouble& param, int n, double variacion)
{
	double valor = 0;
	
	forn(i, n-1)
	{
		param[i] = valor;
		valor += variacion * sqrt(dist(coordenadas[i], coordenadas[i+1]));
	}
	
	/* n siempre es mayor a 0 */
	param[n-1] = valor;
		
	return;
}

/* Deja el parametro t entre 0 y 1. */
void enRango(vdouble& param, int n)
{
	double fin = param[n-1];
	
	fornx(i, 1, n)
		param[i] /= fin;
	
	return;
}

void param_puntos(vPunto& coordenadas, vdouble& param, vPunto& puntos_x, vPunto& puntos_y, int n)
{
	forn(i, n)
	{
		puntos_x[i].x = param[i];
		puntos_x[i].y = coordenadas[i].x;
		
		puntos_y[i].x = param[i];
		puntos_y[i].y = coordenadas[i].y;
	}
	
	return;
}




void muestreo_Uniforme(vdouble& muestreo, int m)
{
	double cantidad = 1.0/(m-1);
	double valor = 0;	
	
	forn(i, m)
	{
		muestreo[i] = valor;
		valor += cantidad;
	}
	
	return;
}
\end{lstlisting}
\end{subsubsection}

\begin{subsubsection}{spline}
\begin{lstlisting}
#include "spline.h"
#include "evaluar.h"

void interpolar(vPunto& puntos, vSpline& spline, int n)
{
	/* n puntos, n-1 splines. */ 
	
	vdouble diag (n-2);
	vdouble diag_otra (n-1);
	/* Representa la diagonal de arriba y de abajo, que son iguales. */	
	
	vdouble res (n-2);
	vdouble c (n);
	
	vdouble h (n-1);
	vdouble a (n-1);
		 
	/* Llenar los vectores. */
	forn(i, n-1)
	{
		h[i] = puntos[i+1].x-puntos[i].x;		
	}
	
	forn(i, n-2)
	{
		diag[i] = 2*(h[i]+h[i+1]);
		diag_otra[i] = h[i];
	}
	
	diag_otra[0] = 0;
	diag_otra[n-2] = 0;
		
	forn(i, n-1)
	{
		/* a[i] = f(x_i+1) - f(x_i). */
		a[i] = puntos[i+1].y-puntos[i].y;	
	}	
	
	forn(i, n-2)
	{
		res[i] = 3.0/h[i+1] * a[i+1] - 3.0/h[i] * a[i];		
	}	

	/* Resolver el sistema con la matriz tridiagonal representada 
	 * por los vectores diag, diag_otra. 
	 * Se usa el algoritmo de matriz tridiagonal.*/
	
	vdouble coef_c (n-2);
	vdouble coef_d (n-2);
	
	coef_c[0] = diag_otra[1]/diag[0];
	coef_d[0] = res[0]/diag[0];
	
	fornx(i, 1, n-2)
	{
		coef_c[i] = diag_otra[i+1]/ (diag[i] - coef_c[i-1] * diag_otra[i]);
		coef_d[i] = (res[i] - coef_d[i-1] * diag_otra[i])/ 	(diag[i] - coef_c[i-1]*diag_otra[i]);
	}
	
	vdouble x (n-2);
	
	x[n-3] = coef_d[n-3];
	fornd(i, n-3)
	{
		x[i] = coef_d[i] - coef_c[i] * x[i+1];
	}	
	
	c[0] = 0;
	c[n-1] = 0;
	
	forn(i, n-2)
	{
		c[i+1] = x[i];
	}
	
	/* Obtener el resto de los coeficientes en funcion de los c_i. */	
	vdouble b (n-1);
	vdouble d (n-1);
	
	forn(i, n-1)
	{
		d[i] = (c[i+1] - c[i]) / (3.0 * h[i]);
		b[i] = a[i]/h[i] - h[i]/3.0 * (2*c[i] + c[i+1]);
	}
	
	/* Crear el spline con los valores ya calculados. */
	forn(i, n-1)
	{
		spline[i].x_i = puntos[i].x;
		
		spline[i].min = puntos[i].x;
		spline[i].max = puntos[i+1].x;
		
		/* S(xj) = f (xj)*/
		
		spline[i].a = puntos[i].y;
		spline[i].b = b[i];
		spline[i].c = c[i];
		spline[i].d = d[i];
	}
	
	return;
}

/* (x-S_x(t))^2 + (y - S_y(t))^2 */
double distancia(Punto& click, vSpline& spline_x, vSpline& spline_y, double punto, int n)
{
	double dist_x = (click.x - evaluar(punto, spline_x, n));
	double dist_y = (click.y - evaluar(punto, spline_y, n));
		
	return (dist_x * dist_x + dist_y * dist_y);
}

/* (S_x(t)-x)((S_x)'(t)) + (S_y(t)-y)((S_x)'(t)). */
double derivada(Punto& click, vSpline& spline_x, vSpline& spline_y, double punto, int n)
{
	double res = (click.x - evaluar(punto, spline_x, n))*(-evaluar_derivada_primera(punto, spline_x, n));
	res += (click.y - evaluar(punto, spline_y, n))*(-evaluar_derivada_primera(punto, spline_y, n));
	
	return res;
}

/* (-(S_x)'(t))^2 + (x-S_x(t))*(-(S_x)''(t)) + (-(S_y)'(t))^2 + (y-S_y(t))*(-(S_y)''(t)). */
double derivada_seg(Punto& click, vSpline& spline_x, vSpline& spline_y, double punto, int n)
{
	double res = ( ((-evaluar_derivada_primera(punto, spline_x, n)) * (-evaluar_derivada_primera(punto, spline_x, n))) + ((evaluar(punto, spline_x, n) - click.x) * (-evaluar_derivada_segunda(punto, spline_x, n))) + ((-evaluar_derivada_primera(punto, spline_y, n)) * (-evaluar_derivada_primera(punto, spline_y, n))) + ((evaluar(punto, spline_y, n) - click.y) * (-evaluar_derivada_segunda(punto, spline_y, n))) );
	
	return res;
}

double mas_Proximo(long int cant_puntos, vSpline& spline_x, vSpline& spline_y, Punto& click, Punto& proximo, int n)
{
	/* f(t) = (x-S_x(t))^2 + (y - S_y(t))^2
	 * f'(t) = 2(x-S_x(t))(-(S_x)'(t)) + 2(y - S_y(t))(-(S_x)'(t)) 
	 * f'(t) = 0 => (S_x(t)-x)((S_x)'(t)) + (S_y(t)-y)((S_x)'(t)) = 0 */	 
	
	vdouble t_partes (cant_puntos);
	muestreo_Uniforme(t_partes, cant_puntos);
	vdouble mejores_t (cant_puntos);
	
	double t_biseccion = multi_biseccion(cant_puntos, t_partes, mejores_t, spline_x, spline_y, click, n);	
	double t_regula_falsi = multi_regula_falsi(cant_puntos, t_partes, mejores_t, spline_x, spline_y, click, n);
	cout << endl << "t de biseccion: " << t_biseccion << " t de regula falsi: " << t_regula_falsi << endl;
	
	double dist_biseccion = fabs(distancia(click, spline_x, spline_y, t_biseccion, n));
	double dist_regula = fabs(distancia(click, spline_x, spline_y, t_regula_falsi, n));
	
	double t;
	/* Veo cual es realmente el mas chico */
	if (dist_biseccion <= dist_regula)
	{
		t = t_biseccion;
	}
	else
	{
		t = t_regula_falsi;
	}
	
	t = newton(t, click, spline_x, spline_y, n);
	
	cout << "t de newton: " << t << endl;
	
	proximo.x = evaluar(t, spline_x, n);
	proximo.y = evaluar(t, spline_y, n);
	
	return t;
}

double biseccion(double a, double b, Punto& click, vSpline& spline_x, vSpline& spline_y, int n)
{
	/* Si la funcion no cambia de signo es razonable suponer que el cero no esta en ese intervalo, 
	 * si el intervalo es lo suficientemente chico. */
	if( (derivada(click, spline_x, spline_y, a, n)) * (derivada(click, spline_x, spline_y, b, n)) > 0)
	{
		if ( fabs(derivada(click, spline_x, spline_y, a, n)) < fabs(derivada(click, spline_x, spline_y, b, n)))
			return a;
		else
			return b;
	}
	
	/* Me fijo en los bordes. */
	if (fabs(derivada(click, spline_x, spline_y, a, n)) < epsilon)
		return a;
		
	if (fabs(derivada(click, spline_x, spline_y, b, n)) < epsilon)
		return b;
	
	double c = (a+b)/2;
	int cant = 0;
	double aux = (derivada(click, spline_x, spline_y, c, n));
	
	while( fabs(aux) > epsilon && cant < MAX )
	{
		if ((derivada(click, spline_x, spline_y, a, n) * aux < 0))
			b = c;
		else
			a = c;
		
		c = (a+b)/2;
		aux = (derivada(click, spline_x, spline_y, c, n));
		cant++;
	}
		
	return c;
}

double multi_biseccion(long int cant_puntos, vdouble& t_partes, vdouble& mejores_t, vSpline& spline_x, vSpline& spline_y, Punto& click, int n)
{ 
	double t = 0;
	
	forn(i, cant_puntos-1) 
	{
		t = biseccion(t_partes[i], t_partes[i+1], click, spline_x, spline_y, n); 
		mejores_t[i] = t;
	}
	
	double dist_min = fabs(distancia(click, spline_x, spline_y, mejores_t[0], n));
	double dist_aux;
	
	/* Veo cual es realmente el mas chico */
	forn(i, cant_puntos-1)
	{
		dist_aux = fabs(distancia(click, spline_x, spline_y, mejores_t[i], n));
		
		if (dist_aux < dist_min)
		{
			dist_min = dist_aux;
			t = mejores_t[i];		
		}
	}
	
	return t;
}

/* Metodo de Newton */ 
double newton(double t_0, Punto& click, vSpline& spline_x, vSpline& spline_y, int n)
{	
	//  x_n+1 = x_n - (f(x_0) / f'(x_0))
	
	double t_res = t_0;
	int cant = 0;
	double aux = derivada(click, spline_x, spline_y, t_res, n);
	
	while( fabs(aux) > epsilon && cant < MAX )
	{		
		/*calculo f'(t)*/
		double aux_deriv = derivada_seg(click, spline_x, spline_y, t_res, n);
		
		/*calculo el x_n+1*/
		t_res = t_res - ( aux / aux_deriv);
		cant++;
		
		/*calculo f(t) */
		aux = derivada(click, spline_x, spline_y, t_res, n);
	};
	
	return  t_res;      
}

/* Metodo de Regula Falsi. */
double regula_falsi(double t_0, double t_1, Punto& click, vSpline& spline_x, vSpline& spline_y, int n)
{
	double t_n_1 = t_0;
	double t_n = t_1;
	double t_res = 0;	
	
	double aux1 = derivada(click, spline_x, spline_y, t_n, n);
	double aux2 = derivada(click, spline_x, spline_y, t_n_1, n);
	 
	if( (aux1 * aux2) > 0)
	{
		if ( fabs(aux1) < fabs(aux2) )
			return t_n;
		else
			return t_n_1;
	}	
	
	/* Me fijo en los bordes. */
	if (fabs(aux1) < epsilon)
		return t_n;
		
	if (fabs(aux2) < epsilon)
		return t_n_1;
	
	double diff = fabs(t_n - t_n_1);
	int cant = 0;
	t_res = t_n - ((aux1*diff) / (aux1-aux2));
	double aux3 = (derivada(click, spline_x, spline_y, t_res, n));
	
	while( fabs(aux3) > epsilon && cant < MAX )
	{
		if(aux3*aux1 < 0)
		{
			t_n = t_res; 
			t_n_1 = t_n;

		}
		else if(aux3 * aux2 < 0)
		{
			t_n = t_res; 
		}
		
		aux1 = (derivada(click, spline_x, spline_y, t_n, n));		
		aux2 = (derivada(click, spline_x, spline_y, t_n_1, n));
		
		diff = fabs(t_n - t_n_1);
		
		/* Calculo x_{n+1}. */
		t_res = t_n - ((aux1*diff) / (aux1-aux2));
		aux3 = (derivada(click, spline_x, spline_y, t_res, n));	
		
		cant++;	
	}
	
	return t_res;		 
}

double multi_regula_falsi(long int cant_puntos, vdouble& t_partes, vdouble& mejores_t, vSpline& spline_x, vSpline& spline_y, Punto& click, int n)
{ 
	double t = 0;
	
	forn(i, cant_puntos-1) 
	{
		t = regula_falsi(t_partes[i], t_partes[i+1], click, spline_x, spline_y, n); 
		mejores_t[i] = t;
	}
	
	double dist_min = fabs(distancia(click, spline_x, spline_y, mejores_t[0], n));
	double dist_aux;
	
	/* Veo cual es realmente el mas chico */
	forn(i, cant_puntos-1)
	{
		dist_aux = fabs(distancia(click, spline_x, spline_y, mejores_t[i], n));
		
		if (dist_aux < dist_min)
		{
			dist_min = dist_aux;
			t = mejores_t[i];		
		}
	}
	
	return t;
}



int ubicacion(Punto& nuevo, vPunto& puntos)
{
	int pos = 0;
	
	while(puntos[pos].x < nuevo.x)
		pos++;
		
	return pos;
}

void agregar_punto(Punto& nuevo, int posicion, vPunto& puntos, vPunto& nuevosPuntos)
{
	/* Paso los puntos a una lista. */
	lPunto puntosAux;
	std::copy (puntos.begin (), puntos.end (), std::back_inserter (puntosAux)); 
	
	/* Inserto en nuevo punto donde corresponde */
	lPunto::iterator it = puntosAux.begin();
	forn(i, posicion)
		it++;
	puntosAux.insert(it, nuevo);
	
	/* Pasar los puntos a un vector nuevamente. */
	std::copy (puntosAux.begin (), puntosAux.end (), std::back_inserter (nuevosPuntos));
	
	return;
}

void graficar(long int cant_puntos, string nombre_fig, string caption, vSpline& spline_x, vSpline spline_y, int n)
{
	vdouble puntos_evaluacion (cant_puntos);
	muestreo_Uniforme(puntos_evaluacion, cant_puntos);
	vPunto puntos_grafico (cant_puntos);
	evaluar(puntos_evaluacion, puntos_grafico, spline_x, spline_y, n, cant_puntos);
	
	ofstream inputGrafico ("input_grafico.m");
	
	stringstream aux_x;
	stringstream aux_y;
	aux_x << "x = [";
	aux_y << "y = [";		
			
	forn(i, cant_puntos)
	{
		aux_x << puntos_grafico[i].x << " ";			
		aux_y << puntos_grafico[i].y << " ";			
	}
	
	aux_x << "];";
	aux_y << "];";
	
	/* Formato matlab. */
	string x, y;
	x = aux_x.str();
	y = aux_y.str();
	
	inputGrafico << "close all" << endl;
	inputGrafico << x << endl;
	inputGrafico << y << endl;
	inputGrafico << "plot(x,y)" << endl;
	inputGrafico << "xlabel('x')" << endl;
	inputGrafico << "ylabel('y')" << endl;
	inputGrafico << "title('" << caption << "')" << endl;
	inputGrafico << "print -deps " << nombre_fig << ".eps" << endl;
	
	//~ system("matlab < input_grafico.m");
	system("octave --quiet < input_grafico.m");
	
	stringstream aux;
	aux << "mv -f " << nombre_fig << ".eps ./graficos";	
	system(aux.str().c_str());
	
	return;
}
\end{lstlisting}
\end{subsubsection}


\begin{subsubsection}{verifar}
\begin{lstlisting}
#include "evaluar.h"

bool igualdad_valida(double a, double b)
{
	return (fabs(a-b) <= epsilon);
}

bool verificar(vPunto& puntos, vSpline& spline, int n)
{	
	forn(i, n-1)
	{
		/* Que interpole bien. */
		if (!igualdad_valida(spline[i].a, puntos[i].y))
		{
			cerr << "No interpola bien en " << puntos[i].x << " " << puntos[i].y << " i: " << i << endl;
			return false;
		}
	}
		
	forn(i, n-2)
	{
		/* Valores iguales de ambos lados. */
		double punto = spline[i+1].x_i;		
		
		if (!igualdad_valida(evaluar_fijo(punto, spline[i]), evaluar_fijo(punto, spline[i+1])))
		{
			cerr << "Los valores no son iguales de ambos lados " << punto << " " << " i: " << i << endl;
			cerr << evaluar_fijo(punto, spline[i]) << " " << evaluar_fijo(punto, spline[i+1]) << endl;
			return false;
		}
		
		if (!igualdad_valida(evaluar_derivada_primera(punto, spline[i]), evaluar_derivada_primera(punto, spline[i+1])))
		{
			cerr << "Las derivadas primeras no son iguales de ambos lados " << punto << " " << " i: " << i << endl;
			cerr << evaluar_derivada_primera(punto, spline[i]) << " " << evaluar_derivada_primera(punto, spline[i+1]) << endl;
			return false;
		}
			
		if (!igualdad_valida(evaluar_derivada_segunda(punto, spline[i]), evaluar_derivada_segunda(punto, spline[i+1])))
		{
			cerr << "Las derivadas segundas no son iguales de ambos lados " << punto << " " << " i: " << i << endl;
			cerr << evaluar_derivada_segunda(punto, spline[i]) << " " << evaluar_derivada_segunda(punto, spline[i+1]) << endl;
			return false;
		}
	}
	
	/* Spline natural. */
	if (!igualdad_valida(spline[0].c, 0))
	{
		cerr << "No cumple la condicion de borde en x_o " << spline[0].c << endl;
		return false;
	}
			
	if (!igualdad_valida(evaluar_derivada_segunda(puntos[n-1].x, spline[n-2]), 0))
	{
		cerr << "No cumple la condicion de borde en x_n " << evaluar_derivada_segunda(spline[n-1].max, spline[n-1]) << endl;
		return false;
	}
		
	return true;
}


double evaluar_fijo(double punto, Spline& spline)
{
	/* Recibe el S_i correspondiente y evalua en punto. */
	double res = 0;	
	
	/* d(x-x_i)^3 + c(x-x_i)^2 + b(x-x_i) + a */
	double aux = punto - spline.x_i;	
	res = spline.d * aux * aux * aux + spline.c * aux * aux + spline.b * aux + spline.a;
	
	return res;
}

double evaluar_derivada_primera(double punto, Spline& spline)
{
	/* Recibe el S_i correspondiente y evalua en punto. */
	double res = 0;	
	
	/* 3d(x-x_i)^2 + 2c(x-x_i) + b */
	double aux = punto - spline.x_i;
	res = spline.d * 3 * aux * aux + spline.c * 2 * aux + spline.b;
	
	return res;
}

double evaluar_derivada_segunda(double punto, Spline& spline)
{
	/* Recibe el S_i correspondiente y evalua en punto. */
	double res = 0;	
	
	/* 6d(x-x_i) + 2c */
	double aux = punto - spline.x_i;			
	res = spline.d * 6 * aux + spline.c * 2;
	
	return res;
}

double evaluar(double s, vSpline& spline, int n)
{
	double res = 0;
	
	/* Encontrar el S_i correspondiente y evaluar. */
	forn(i, n-1)
	{
		if (spline[i].min <= s && s <= spline[i].max)
		{
			/* a(x-x_i)^3 + b(x-x_i)^2 + c(x-x_i) + d */
			double aux = s - spline[i].x_i;			
			res = spline[i].d * aux * aux * aux + spline[i].c * aux * aux + spline[i].b * aux + spline[i].a;
			break;
		}
	}
	
	return res;
}

double evaluar_derivada_primera(double s, vSpline& spline, int n)
{
	double res = 0;
	
	/* Encontrar el S_i correspondiente y evaluar. */
	forn(i, n-1)
	{
		if (spline[i].min <= s && s <= spline[i].max)
		{
			res = evaluar_derivada_primera(s, spline[i]);
			break;
		}
	}
	
	return res;
}




double evaluar_derivada_segunda(double s, vSpline& spline, int n)
{
	double res = 0;
	
	/* Encontrar el S_i correspondiente y evaluar. */
	forn(i, n-1)
	{
		if (spline[i].min <= s && s <= spline[i].max)
		{
			res = evaluar_derivada_segunda(s, spline[i]);
			break;
		}
	}
	
	return res;
}


void evaluar(vdouble& muestreo, vPunto& muestreados, vSpline& spline_x, vSpline& spline_y, int n, int m)
{
	double X, Y, s;
	
	forn(i, m)
	{		
		s = muestreo[i];		
		X = evaluar(s, spline_x, n);
		Y = evaluar(s, spline_y, n);				
		
		muestreados[i].x = X;
		muestreados[i].y = Y;		
	}
	
	return;
}
\end{lstlisting}
\end{subsubsection}


\end{subsection}

%~ \newpage
%~ \begin{subsection}{Imágenes}

%~ \begin{figure}[H] \centering
	%~ \subfloat{\includegraphics[scale=0.3]{./otros/Ramanujan.png}}	
	%~ \caption{Ramanujan}
%~ \end{figure}

%~ \end{subsection}

\end{section}
